Analyze ZFP503+ cells from Macosko et al. single cell mouse retina data

Introduction

David, Please see the link below for the paper on Single-Cell retinal sequencing paper. We can ask the following questions. Whether Nlz2 present one or more clusters? Whether Nlz2 present in Amacrine or ganglion cell clusters? If present, what other transcriptional factors present in those clusters?

http://www.sciencedirect.com/science/article/pii/S0092867415005498

There is another paper on Single-cell sequencing done on Bipolar cell. See the link below. http://www.sciencedirect.com/science/article/pii/S0092867416310078

Thanks for your help.

Elan

Visualization of single cell expression, clustered by t-SNE

Seurat was used for analysis.

Briefly, the raw counts (P14 mouse retina) data from the Macosko and McCarroll paper were loaded into Seurat and, for visualization purposes, only the 12934 cells with 900 or more genes expressed in them were used for initial clustering and visualization. This cut-off was used as the majority of the cells have very little detected gene expression and the clustering was not performing well.

My clustering results

library(cowplot)
library(Seurat)
library(tidyverse)

Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
package ‘tibble’ was built under R version 3.4.1package ‘tidyr’ was built under R version 3.4.1package ‘purrr’ was built under R version 3.4.1package ‘dplyr’ was built under R version 3.4.1Conflicts with tidy packages -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
accumulate(): purrr, foreach
combine():    dplyr, Biobase, BiocGenerics
expand():     tidyr, Matrix
filter():     dplyr, stats
ggsave():     ggplot2, cowplot
lag():        dplyr, stats
lift():       purrr, caret
Position():   ggplot2, BiocGenerics, base
when():       purrr, foreach

library(here)
# load processed data from scripts/process_macosko.R
load(here('data/retina_seurat_subSet.Rdata'))
# Choose gene outliers on mean-variability plot
#retina <- FindVariableGenes(object = retina, x.low.cutoff = 0, y.cutoff = 2)
#length(x = retina@var.genes)
# elbow at 30
#PCElbowPlot(object = retina, num.pc = 40)
#PrintPCA(object = retina, pcs.print = 1:36)
#PCHeatmap(object = retina, pc.use = 1:12,100)
#PCHeatmap(object = retina, pc.use = 13:24,100)
#PCHeatmap(object = retina, pc.use = 25:36,100)
#VizPCA(object = retina, pcs.use = 1:2)
TSNEPlot(object = retina, do.label = TRUE)

As Macosko et al. have done the difficult work of classifying the cell types already, I will be using their cell classifications from now on.

Macosko Cluster and Cell Type Table/Key

To make sure that my clustering approximately matches the Macosko assignments, I have replotted the t-SNE plot from above, this time with the Macosko cell group assignments. They look pretty good, so it seems reasonable to use them.

Cell clustering with Macosko et al. assignments

# remove cells that Macosko did not classify
notNA <- retina@meta.data[!is.na(retina@meta.data$Macosko_Clusters),] %>% rownames()
retina_s <- SubsetData(retina, cells.use = notNA)
#print with macosko clustering
TSNEPlot(object = retina_s, do.label = TRUE, group.by = 'Macosko_Clusters')

Plot retina marker genes and ZFP503.

I can add any other genes you would like to see.

Two points:

1. Not that many ZFP503 cells
2. Many of them are in cluster 9 (top), which is an Amacrine population

# plot individual genes
FeaturePlot(object = retina, features.plot = c("PAX6", "VSX2", "OPN1MW", "RLBP1","NRL","SLC17A6","CRX","ZFP503"), 
            pt.size = 0.1,
            cols.use = c("gray", 
                         "blue"))

Number of cells with detectable ZFP503 Expression.

828 cells out of 44098 have detectable ZFP503 expression. Not many. This doesn’t seem to be a very highly expressed gene. This is corroborated by looking at ZNF503 expression in eyeIntegration

table(retina_superset@data['ZFP503',] > 0.1)

FALSE  TRUE 
43270   828 

Expression of key retinal markers in ZFP503 positive cells

Macosko Clustering cell type assignments for ZFP503

Amacrine enrichment. Also in the rods, but at a lower proportion than you would expect by chance.

Differential Expression

Find differentially expressed genes in ZFP503 positive cells relative to non-ZFP503 expressing cells. pct.1 is the percentage of ZFP503+ cells that the gene is expressed in. pct.2 are the ZFP503- cells. The p_val is calculated off of the avg_diff stat, which measures the strength of expression. So you will see genes that don’t have a large difference in percentage of cells expression (pct.1 and pct.1), but have a very low p_val. A positive avg_diff means that the gene is more highly expressed in the ZFP503+ cells.

SessionInfo

sessionInfo()

R version 3.4.0 (2017-04-21)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2        dplyr_0.7.3         purrr_0.2.3         readr_1.1.1         tidyr_0.7.1         tibble_1.3.4        tidyverse_1.1.1     Seurat_2.0.1        Biobase_2.36.2      BiocGenerics_0.22.0 Matrix_1.2-11      
[12] cowplot_0.8.0       doMC_1.3.4          iterators_1.0.8     foreach_1.4.3       caret_6.0-77        ggplot2_2.2.1       lattice_0.20-35     here_0.1           

loaded via a namespace (and not attached):
  [1] readxl_1.0.0         backports_1.1.0      Hmisc_4.0-3          VGAM_1.0-4           NMF_0.20.6           sn_1.5-0             plyr_1.8.4           igraph_1.1.2         lazyeval_0.2.0       splines_3.4.0       
 [11] gridBase_0.4-7       digest_0.6.12        htmltools_0.3.6      lars_1.2             gdata_2.18.0         magrittr_1.5         checkmate_1.8.3      cluster_2.0.6        doParallel_1.0.10    mixtools_1.1.0      
 [21] ROCR_1.0-7           recipes_0.1.0        modelr_0.1.1         gower_0.1.2          dimRed_0.1.0         R.utils_2.5.0        colorspace_1.3-2     rvest_0.3.2          haven_1.1.0          jsonlite_1.5        
 [31] lme4_1.1-13          bindr_0.1            survival_2.41-3      ape_4.1              glue_1.1.1           DRR_0.0.2            registry_0.3         gtable_0.2.0         ipred_0.9-6          MatrixModels_0.4-1  
 [41] kernlab_0.9-25       ddalpha_1.2.1        prabclus_2.2-6       DEoptimR_1.0-8       scales_0.5.0         mvtnorm_1.0-6        rngtools_1.2.4       Rcpp_0.12.12         dtw_1.18-1           xtable_1.8-2        
 [51] htmlTable_1.9        tclust_1.3-1         foreign_0.8-69       proxy_0.4-17         mclust_5.3           SDMTools_1.1-221     Formula_1.2-2        DT_0.2               tsne_0.1-3           stats4_3.4.0        
 [61] lava_1.5             prodlim_1.6.1        httr_1.3.1           htmlwidgets_0.9      FNN_1.1              gplots_3.0.1         RColorBrewer_1.1-2   fpc_2.1-10           acepack_1.4.1        modeltools_0.2-21   
 [71] ica_1.0-1            pkgconfig_2.0.1      R.methodsS3_1.7.1    flexmix_2.3-14       nnet_7.3-12          labeling_0.3         rlang_0.1.2          reshape2_1.4.2       cellranger_1.1.0     munsell_0.4.3       
 [81] tools_3.4.0          ranger_0.8.0         broom_0.4.2          evaluate_0.10.1      stringr_1.2.0        yaml_2.1.14          ModelMetrics_1.1.0   knitr_1.17           robustbase_0.92-7    caTools_1.17.1      
 [91] pbapply_1.3-3        nlme_3.1-131         R.oo_1.21.0          RcppRoll_0.2.2       xml2_1.1.1           compiler_3.4.0       e1071_1.6-8          ggjoy_0.3.0          stringi_1.1.5        forcats_0.2.0       
[101] trimcluster_0.1-2    psych_1.7.5          nloptr_1.0.4         diffusionMap_1.1-0   data.table_1.10.4    bitops_1.0-6         irlba_2.2.1          R6_2.2.2             latticeExtra_0.6-28  KernSmooth_2.23-15  
[111] gridExtra_2.2.1      codetools_0.2-15     MASS_7.3-47          gtools_3.5.0         assertthat_0.2.0     CVST_0.2-1           pkgmaker_0.22        rprojroot_1.2        withr_2.0.0          mnormt_1.5-5        
[121] hms_0.3              diptest_0.75-7       grid_3.4.0           rpart_4.1-11         timeDate_3012.100    minqa_1.2.4          class_7.3-14         rmarkdown_1.6        segmented_0.5-2.1    Rtsne_0.13          
[131] numDeriv_2016.8-1    scatterplot3d_0.3-40 lubridate_1.6.0      base64enc_0.1-3